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Abstract 

To assess the role of single-chain elasticity, non-affine strain helds and pre¬ 
stressed reference states we present and discuss the results of numerical and 
analytical analyses of modihed 8-chain Arruda-Boyce model for cross-linked 
polymer networks. This class of models has proved highly successful in mod¬ 
eling the hnite-strain response of flexible rubbers. We extend it to include 
the effects of spatial disorder and the associated non-affinity, and use it to 
assess the validity of replacing the constituent chain’s nonlinear elastic re¬ 
sponse with equivalent linear, Hookean springs. Surprisingly, we hnd that 
even in the regime of linear response, the full polymer model gives very dif¬ 
ferent results from its linearized counterpart, even though none of the chains 
are stretched beyond their linear regime. We demonstrate that this effect is 
due to the fact that the polymer models are under considerable pre-stress in 
their ground state. We show that pre-stress strongly suppresses non-affinity 
in these unit cell models, resulting in a marked stiffening of the bulk response. 
The effects of pre-stress we discuss may explain why fully affine mechanical 
models, in many cases, predict the bulk mechanical response of disordered 
stiff polymer networks so well. 
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1. Introduction and summary 

The mechanical response of soft polymeric materials such as rubbers, plas¬ 
tics, gels and hlamentous biomaterials is crucial to their technological and 
biological functioning. While it is obvious that the macroscopic response 
of these materials is encoded in the behavior of its constituent chains and 
the manner in which these are connected and arranged in space, predicting 
the mechanical response based on such microstructural parameters remains 
a challenge. Generally, the hlaments themselves display strongly nonlinear 
response to deformations, the spatial arrangement is highly disordered, the 
strain helds and distributions may be highly heterogeneous and hnally, pre¬ 
dictive models far beyond linear response are required to capture typical 
operational conditions for biological and synthetic materials alike. 

Literature abounds with different approaches put forward to tackle one, or 
several, of these challenges: at the microscopic end of the spectrum, molecular 
dynamics and Monte Carlo models allow access to micromechanical response 
using accurate representations of chemical or supramolecular polymer ma¬ 
terial constituents, but are restricted both in the largest sizes accessible, as 
well as in the timescales. At the largest scales, phenomenological constitutive 
models are often able to capture nonlinear bulk response well for specihc ma¬ 
terials, but lack a direct link to microscopic properties and architecture and 
therefore must, at some point, invoke effective parameters determined from 
hts to data. In between the length scales, so-called rubber elasticity mod¬ 
els [D12113] coarse grain the molecular response into effective ideal irnsiE], 
wormlike [7] or differently elastic chain response and, in order to connect 
to macroscopic response, apply averaging or homogenization techniques to 
compute the macroscopic response. Each of these approaches has proven ex¬ 
tremely useful, and each continues to be improved. From a practical point of 
view, provided one selects the appropriate approach for a given system and 
set of mechanical questions, the available spectrum of techniques is eminently 
capable of accurate, even predictive, analysis. From a fundamental point of 
view, however, the situation remains unsatisfactory: there is to date no single 
model that allows predictive modeling of macroscopic mechanical response 
using as input the single-chain response and the microstructural architecture 
of the material. Nor is it clear, in fact, what processes and mechanisms such 
a model should include. Nonlinearities enter at all levels: from the polymers 
themselves, from spatial and elastic heterogeneity, from composition and ge¬ 
ometry. With this paper, we assess the impact of several of these factors - 
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specifically, the elastic nonlinearity of individual polymers, non-affine strain 
fields, and pre-stresses - by introducing them, in controlled fashion, into an 
established framework. 

Our point of departure is a particular class of models: 8-chain or Arruda- 
Boyce (AB) models. We are aware that numerous improved and more in¬ 
tricate versions have been proposed il El El 0 EDI EH 112], but all share the 
ability to determine the constitutive relation directly from single-polymer re¬ 
sponse, for a specihc microstructural arrangement of chains. The AB model 
shares two important features with classical (Flory) rubber elasticity (CRE) 
[T]: It is (macroscopically) isotropic and directly rooted in the response of 
single, freely jointed chains (FJCs). One important advantage of the AB 
model is that it explicitly considers the spatial arrangement of these chains, 
rather than apply to them a statistical averaging as CRE does (its response 
is computed by averaging over all chain orientations, weighted according to 
the FJC’s Gaussian radial distribution function). In contrast to the CRE 
model, the AB is a unit-cell model: It comprises eight identical freely jointed 
chains arranged in symmetric fashion within a cube of linear dimension oq. 
This basic motif readily permits the introduction of spatial disorder, and 
numerous more intricate, improved variants of the model exist |H El E] • hi 
the following, we study the AB model for different types of polymer elastic 
behaviour, with the original - ordered - geometry (chains linked in the center 
of the unit cell) and with positional disorder (chains linked off-center). Our 
work shows, that there is a direct connection between single-chain elasticity, 
disorder, non-affinity and pre-stress that has important consequences for the 
predicted moduli and the response. 

Our paper is organized as follows: in Sec. |^we review the Arruda-Boyce 
model. In Sec. [^we present and discuss our disordered, non-affine extension 
of it. We demonstrate that the modihed model is capable of capturing ex¬ 
perimental data. In Sec. |^we compare the stress-strain response in various 
settings to reveal the importance of pre-stressed states. In Sec. [^we present 
a measure of non-affinity and compare our Endings in both the pre-stressed 
and stress free ensembles. In Sec. [^we present our conclusions and add some 
remarks regarding the relevance of pre-stresses in a more general setting and 
how they afiect the degree of affinity of deformations. 
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2. Freely Jointed Chain mechanics and the Arruda-Boyce model 


In the standard Arruda-Boyce 8-chain model [1] a polymer network is 
represented as a collection of identical unit cells. Each of these cubic unit 
cells has a linear dimension of oq in the undeformed state, and contains 8 
polymer chains that connect the cube’s vertices to its center. Thus, the 8 
chains are identical in initial length £o: 

£o = ^V3ao. ( 1 ) 

Each of these chains is modeled as a Freely Jointed Chain (FJC): the well- 
known model for flexible (Gaussian) polymers [21 HI US]- The FJC represents 
a long polymer with a contour length Got as a sequence of N connected, but 
freely hinged, Kuhn segments each of which has a length b (such that Nb = 
Got), where b is chosen several times larger than the persistence length of the 
polymer to ensure that no orientational correlations between the segments 
remain. 

The FJC has no internal energy - each of its conformations may be trans¬ 
formed into any of the other conformations at no energetic cost. Its Gibbs 
free energy, in the presence of an external applied force / in the ^direction, 
is therefore 

g = H-TS = -fi, - TS , ( 2 ) 

where H is the enthalpy, T is the absolute temperature, S is the entropy and 
iz is the end-to-end length of the polymer in the direction of the force, i.e. it 
is the 5 component of the full end-to-end vector i which - expressed in terms 
of the unit bond vectors ti - is computed as 

N 

i=bY,k. ( 3 ) 

i=l 


The partition function in the Gibbs ensemble may be computed as 

N 


Z{N, f,T)= /dG ■ ■ ■ /dfjv exp 


ksT 


i=l 


i ■ Z 


(4) 


In three dimensions, the partition function is readily computed and from it, 
an analytical expression for the Gibbs free energy is obtained 


G{N,f,T) 


-kBT\ogZ{NJ,T) 

Att smh.{fb/k^T) 


—Nk^T log 


fb/k^T 


(5) 
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We may perform a Legendre transformation on Q{N,f,T) to obtain from 
this the Helmholtz free energy T = Q + fiz to find for the free energy of a 
FJC subject to a force / 


J^{N,iz,T) 


—Nk-^T log 


/ 47rsinh(/6/fcBT)\ 

V fb/kBT J 


+ / 4 . 


( 6 ) 


Evaluating now the condition that W be minimal in equilibrium yields the 
mechanical response of the FJC, summarized as the relationship between its 
projected end-to-end length iz and the externally applied force / [2] 


W 

(Qif) 


0 ^ 


Nb 



NbCi 


ksT 


ksT 

fb . 


(7) 


where we define the first order Langevin function: Ci{x) = cothx — x~^. 
The angular brackets, (•), signify ensemble averages. In the following, we 
will denote by i this ensemble average and will therefore drop the brackets 
as well as the subscript z. While the Langevin function does not have an 
analytical inverse, we may use Eq. Q to formally extract the force-extension 
relationship that we will use throughout this paper 


/W 



( 8 ) 


Since we shall generally specify the state of strain, rather than the stress, it 
is convenient to substitute the expression for f{i) in to the formula for the 
Helmholtz free energy Eq. @, to arrive at the expression also used in |1] 


HNJ.T) 



+log 




sink (A Xm)) 


+J^o{N,T) 


(9) 


We may expand this expression for small t, using the Taylor expansion of 
the inverse Langevin function ~ 3x, to find that 


T{N,i,T) 




( 10 ) 
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From which we see that the elasticity, for small extensions, is Hookean with 
an entropic spring constant equal to 


_ 3/cbT 
“ Nh^ • 


( 11 ) 


The equilibrium length is zero for the FJC. 

We now turn to the description of deformation applied to the unit cell. 
We will start with a deformation that comprises three principal stretches Aj, 
and transforms the boundary according to 


/ Ai 0 0 \ 

A = 0 A2 0 , (12) 

VO 0 As / 


which maps each point on the boundary of the unit cell to = A ■ Fb. 
Since the arrangement of chains within the Arruda-Boyce unit cell is fully 
symmetric, we may infer from the deformation of the boundary also the 
displacement of the central point, where the eight chains meet. This is a 
crucial point: by construction, the Arruda-Boyce unit cell is constrained to 
deform affinely. 

This considerable simplification means that specifying {Ai,A 2 ,A 3 } gives 
full knowledge of the deformed lengths of all eight chains. Moreover, be¬ 
cause of their symmetric arrangement in space, in this deformation all chains 
experience identical changes in length, and the new (deformed) projected 
end-to-end length and the original length is given by 

^(Ai, A 2 , A 3 ) = |A ■ .^ol = “^ (Ai + \ 2 ~i~ A 3 ) ^ 5 (13) 

where oq is the linear dimension of the unit cell and £0 is the initial length 
of the chain. In the AB model, the eight chains are modeled as FJC’s and 
the undeformed lengths are identified with the root of the mean-squared 
end-to-end length 

4 = \fWW^) = ( 14 ) 

This identification has important consequences: the undeformed length is 
not the length at zero force, which is zero as is also evident from Eq. 
Thus, the chains in the AB model need to be stretched in order to reach the 
central point, by a tension which may be computed to be 


fps 


kj^ , 



(15) 
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This gives rise to a pre-stress in the resulting model which we will return 
to in the following. With this definition, however, the length of the chains 
may be expressed in terms of the FJC parameters N and b, and the principal 
stretches Aj, as 


^(-^1; -^2, A3) — (A^ + + A3) ^ 


( 16 ) 


By substituting Eq. IT into Eq. we now have an expression for the free 
energy density of deformation expressed in terms of the stretches Aj and the 
FJC parameters. 


/(A,) = pT(N, 4 ViVt (a; + A= + Ai)*'", T). 


(17) 


where p = S/oq is the chain density for the eight-chain AB model in the cube 
of side oq. From this, we may compute the principal (Cauchy) stress tensor 
components 


cr, 


C 


1 

det A 



A. A 

A1A2A3/ 


d\i 


+ p. 


( 18 ) 


Where p is an arbitrary hydrostatic pressure, which may be applied to en¬ 
sure the pre-stresses are compensated for in the reference configuration. For 
incompressible materials, the determinant of A equals one and we hnd, by 
direct differentiation, and after some algebra 


o? = 


5 (pioT) ^^N ( - 


X 


(19) 


The combination pk^T is commonly denoted Cr, the classical rubber elastic 
modulus (although we note that in their original paper, Arruda and Boyce 
denote by CR the full prefactor ^pk^T). In this paper, we will follow con¬ 
ventions and report nominal, rather than Cauchy stresses - these are the 
forces per unit underformed area, which we shall denote by simply a. For an 
incompressible material they may be straightforwardly obtained from by 
correcting for the incurred stretches as 

( 20 ) 
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3. Model definitions and simulation protocol 


Throughout this paper we will analyze and compare nominal stresses for 
various types of homogeneous deformations: 

• uniaxial extension or compression, for which the deformation gradient 
tensor is 

/AO 0 \ 

Ac; = 0 A-3 0 (21) 

Vo 0 A-i / 

• biaxial extension or compression, for which 

/ A-2 0 0 A 

As = 0 AO (22) 

V 0 0 A / 

• simple shear, for which 

/I 0 A\ 

A5 = 0 1 0 (23) 

Vo 0 1 / 


• pure shear (or plane strain compression) for which 


Ap — 


A 

0 

0 


0 

1 

0 


0 

0 

1 

A 


(24) 


• uniform expansion or compression. This deformation mode does not 
conserve volume and we will use to assess the hydrostatic pre-stress 
which is probed by changing the volume. 

/ A 0 0 A 

Ay = 0 A 0 (25) 

V 0 0 A / 


We point out that Eqs. 22 21 are in fact inter-changeable provided that 
A“^ —)■ A or A “2 —>■ a, but for the sake of consistency with previous work we 
will keep both definitions from now on. Our first objective is to use the AB 







Figure 1: Schematic two dimensional representation of a typical Arruda-Boyce setting 
(left) and our model incorporating positional disorder with the off-center linker (right). 


modeling approach to assess the dependence of the predicted mechanical re¬ 
sponse on the type of elasticity chosen for the individual polymers: Hookean 


springs with spring constants as dehned in Eq. 11, and Freely Jointed Chains 
with energy is given by Eq. For the linear elastic Hookean springs, as is 
commonly done in polymer network modeling [H [la [ig HU , we assume a 
rest length equal to the underformed length. The energy of linear elastic 
spring is equal to Taking the linearized spring constant of freely 


jointed chain model for hnite deformations (Eq. 11) we arrive at: 


J^iN I T) K,- ( 

^ ^ 2ViV62 


5t 


(26) 


with 6i now the extension away from the rest length Eq. In the following we 
test whether replacing the nonlinear force-extension of the FJC model with 
its linear, hnite rest-length approximation affects the linear (and nonlinear) 
response of an AB-type network. 

The second objective is to study the effects of allowing for non-afhne 
deformations. In order to do this, we modify the geometry, placing the linker 
connecting the eight chains off-center (this is schematically shown in a 2D 

plot in Fig. [^. 

This introduces important degrees of freedom into the AB model: the 
asymmetric (off-center) placement of the linker allows non-afhne deforma¬ 
tions of the connected chains. That is, if its initial position is ry, its position 
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after the boundaries of the cube have been subjected to A is not A-Vc- Rather, 
the new (deformed) position is determined by minimizing the sum of the elas¬ 
tic energies from all eight chains. The elastic energies thus obtained are then 
averaged among all the disordered unit cells considered. This introduces 
elastic inhomogeneity within the unit cell. When the linker is not centrally 
placed, some chains are shorter than others and therefore stiffer than others: 
As may be surmised from Eq. 10, the spring constant for entropic FJC’s 


scales as N~^, and longer chains are thus less stiff (lower spring constants 
imply less force is required for the same extension) than shorter ones. We 
shall term the hrst effect positional disorder, and the second elastic disorder. 
One hnal effect that is unavoidably incurred is related to the density. Eight 
chains connected with a linker positioned exactly in the center of the unit cell 
requires the least amount of spring contour length, and thus conhgurations 
that are linked elsewhere inside the volume require more spring length per 
unit volume. 

In what follows, we focus on elastic and positional disorder. For the 
polymer models we are considering, the two are inseparable, as the length 
differences that give rise to positional disorder are directly related to stiffness 
differences. In order to disentangle these from density effects, we have de¬ 
vised a computational method that allows one to work in a constant density 
ensemble. For FJCs this has important consequences for mechanical response 
as we will show. 

Previous work has examined the effects of either non-affinity of the bound¬ 
aries (modeling the compliance of the surroundings of the quasi-unit cell) 
|18j . frozen positional disorder (the case where the interior linker point is 
distributed throughout the unit cell, but still constrained to deform affinely) 
[3], and coarse-grained (homogenized) approaches to the average non-affine 
displacements [6]. These works emphasize the large effects of both types of 
disorder on mechanical response, but each is unable to disentangle these ef¬ 
fects from either the constitutive behavior of the chains themselves, or from 
density effects. This is the central novelty of our work: by operating in 
a constant density ensemble, averaging over well-dehned discrete networks, 
we are able to clearly establish the contributions of polymer response and 
non-affinity. 

The constant density ensemble is constructed as follows. We start from 
the observation that each AB unit cell, as dehned in the original paper, 
contains SxN FJC segments, each of length b. In all of our simulations, we 
set Oo to a value of 1 (we will work in non-dimensional units, scaling lengths 
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by Oq from now on, but note that when htting to data the value of CR or, 
alternatively, pk^T, sets all relevant length- and force scales) and = 50. 
The task at hand is to construct all conhgurations that contain exactly 400 
segments, distributed in all possible manners across the eight chains. We 
do this by hrst partitioning the 400 segments into eight subsets 
constructing FJC’s out of each. We then presume that these eight chains are 
connected at some position ry within the cube by solving the force balance 
equation. rc{{Ni}) is thus dehned by the implicit equation 


i=l,8 


fiirc) = 




j=l,8 




-1 


iijrc) 

bN, 


= 0 . 


(27) 


Here, £i{fc) is the unit vector pointing along polymer i. This force balance 
relation consists of three separate force balances along the three axes. Since 
there are three degrees of freedom (the coordinates fy), this equation has a 
unique solution. 

The force balance criterion described by Eq. ^ can be achieved by fol¬ 
lowing any of the two subsequent routes: hrstly, place the linker at a desired 
position fy and retrieve the set W that satishes Eq. ^ and secondly, choose 
a set Ni, under the constant density condition, and compute the correspond¬ 
ing position fy({W}) of the crosslinker. The former route would be the path 
of choice as it would allow to tune the resulting conhguration to any desired 
spatial distribution. This former route is, unfortunately, mathematically un¬ 
feasible, as there are eight unknowns to be retrieved from a system of at 
most four equations: three equations dictating the force balance criteria for 
each Cartesian coordinate and a fourth equation for the constant density 
ensemble. 

We therefore turn to the latter route and start by randomly picking a set 
Ni, solve Eq. 27, store the position of the crosslinker fy({W}) and repeat 


this procedure to yield a spatial distribution similar to the one graphed in 
the upper panel of Fig. Importantly, though all these equilibrated ground 
state conhgurations with off-center linkers Fc({iVj}) have exactly the same 
density, the differ greatly in the amount of pre-stress. In a single chain 
framework, we dehne pre-stress as the amount of tension needed to maintain 
the chain’s ends at specihc points: one end is undoubtedly tethered at one 
vertex of the unit cell and the other end is placed at fc{{Ni}). Each possesses 
internal forces that sum to zero, but as all of the FJCs are stretched away 
from zero length, the tension in each chain is not equal to zero. As we have 
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Figure 2: The ensemble of equal density networks (left panel) that are generated with the 
force balance equation is further subjected to a quasi-random selection (right panel) to 
yield a smooth spatial distribution of the intersection point. 

seen in the previons, the stiffness of chains increases with their length and to 
avoid physically nnrealistic short chains (which contribnte disproportionally 
to the overall response) we place lower-bonnd cntoff of 10 segments on iVj. 

The representation of linker positions obtained in the described manner 
is shown in Fig. with the npper panel plotting the fnll set of Fc({iVj}) that 
meet the force balance criterion and the lower panel with a qnasi-random 
selection. As one can see, the fnll ensemble is highly textured and the den¬ 
sity of points is not constant throughout the unit cube. We infer that this 
textured structure is an immediate result of the nonlinearities present in the 
force-extension curve of individual chains. As this texture is a result of our 
procedure, and does not reflect a physical property of a random sub volume of 
polymer material, we do not sample randomly from this ensemble but rather 
do so correcting for its density variations: we subdivide the unit cube into 
smaller volumes (one thousand cubical boxes, each having an edge length 
equal to 1/10 of the unit cell’s edge length) and pick randomly ten possible 
configurations from within each randomly chosen sub volume, to distribute 
the sampling points evenly across the unit cube. 
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Table 1: The models compared in this paper 


Model 

Linker Placement 

Affine / N onaffine 

Chain Elasticity 

Polymer Arruda-Boyce (PAB) 

central 

affine by symmetry 

FJC 

Hookean Arruda-Boyce (HAB) 

central 

affine by symmetry 

Hookean 

Polymer Non-Affine (PNA) 

distributed 

non-affine 

FJC 

Hookean Non-Affine (HNA) 

distributed 

non-affine 

Hookean 

Polymer Disordered Affine (PDA) 

distributed 

forced affine 

FJC 

Hookean Disordered Affine (HDA) 

distributed 

forced affine 

Hookean 


Once an eqnilibrinm configuration is chosen, it is subjected to the various 
deformations detailed above. During the deformation, we introduce on last 
feature to the procedure: we either allow the linker to displace non affinely, or 
we force it to deform affinely. Thus, we study six distinct variants of Arruda- 
Boyce type unit cell models, evaluated at identical densities but with or 
without nonlinear stretching, with or without positional/elastic disorder and 
with or without the freedom to deform affinely. These models are collected 
Table 1. We now proceed to analyze the differences and similarities among 
them. 

As a hrst check, we compare the ability of the full non-affine unit cell 
model PNA to that of the original Arruda-Boyce model to capture the data 
of Treloar. The results are shown in Fig. which demonstrates that with 
the appropriate parameter choices, (oq = 0.5945, W = 212 and 6 = 0.1) our 
model is equally capable of describing the data, with an equal amount of 
fit parameters but including the microscopic effects of non-affinity. We note 
that the parameters for which this fit was obtained were CR = O.OOMPa, and 
N = 26.5 [1]. While we are generally interested in the generic behavior of 
modified AB models, we do note that as the AB model is one of the few that 
allows one, in principle, to relate macroscopic stresses directly to microscopic 
response these are not simply £t parameters: they do pertain to the actual 
densities, at least those obtained within the context of the model. The fits 
obtained for the Treloar data, for instance, translate into a unit cube with 
a side of about 5nm, and a link length b = 0.83nm. This is very close to 
molecular length scales and it is highly questionable whether FJC’s remain 
credible representations of rubber-like networks at these length scales. 
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A 


Figure 3: Stresses (T[MPa] versus applied strain A in an uni-axial type of deformation. 
Averaged stresses (Q) in a network assembled by modified disordered unit cells plotted 
against experimental data (■) obtained in m and used in the original Arruda-Boyce 
model ID . Both data sets are graphed alongside the analytical curve (black line) predicted 
by the Arruda-Boyce model [3]. The presence of shorter chains in the statistics increases 
stress values for large deformations. 
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4. Comparing the stress-strain response of the six models 


Each of the six models is subjected to the hve different deformations. We 
plot the measured energies, and the measured stresses - evaluated at precisely 
the same chain densities and, for the distributed linker models, averaged over 
the same ensemble of linker positions, in Figures |4]j^ Plotted are the total 
energies 

8 

= (28) 
i=l 


where the energy contribution of each chain separately, is given 

either by Eq. H or Eq. depending on the chain elastic response, and the 
nominal stress (setting p = 0), 


_ 1 dTtot 

V dX 

We consider one deformation type where volume of the unit cell is not pre¬ 
served - the uniform expansion/compression, where the stress needs to be 
corrected. We use of this infinitesimal deformation mode to establish the 
pre-stress. Note the symmetries in the energy-strain curves (upper panel 
of Fig. [^: PDA, PNA and PAB energies are odd around the ground state 
(A = 1), unlike HDA, HNA, HAB that are even around the same state. In the 
lower panel of the same Fig. [^stress-strain curves of the same pre-stressed 
models, PDA, PNA and PAB show finite stresses at no deformation (A = 1). 
A lot of information is presented in rather compact format in these graphs 
but we distill from it three key observations. 

• Observation 1: In uniform extension, FJC models have nonzero stress 
at zero strain while Hookean models do not. 

• Observation 2: For all cases considered, the PNA and PDA models 
give indistinguishable results. For FJC models, disorder matters but 
affinity/non-affinity does not. 

• Observation 3: Both classes considered throughout this paper (FJC 
and Hookean models) show the same type of behavior: introducing 
disorder leads to a more rigid mechanical response when compared to 
Arruda-Boyce settings. 
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Figure 4: Energies (left column) and stresses (right column) versus applied deformation 
A in the homogenous deformations used. Ijigthe left column, we show the energy-strain 
curves for the three polymer (FJC) models: the PDA model (<)), the PNA model (a) and 
the PAB model (*). The inset shows the three Hookean models: the HDA model (v), the 
HNA model (o) and the HAB model (□). The right column shows nominal stress-strain 
curves for the PDA (<)), the PNA (a), the PAB (*), the HDA (v), the HNA (o) and HAB 
models (□). All models are evaluated at 8N = 400, b = 0.1, oq = 1, fesT = 1. 
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Figure 5: Energies (upper panel) and stresses (lower panel) versus applied strain A in 
uniform extension deformation Aj/. In the upper panel, we show the energy-strain 
curves for the three polymer (FJC) models: the PDA model (0), the PNA model (a) 
and the PAB model (*). The upper panel inset shows the three Hookean models: the 
HDA model (v), the HNA model (o) and the HAB model (□). The lower panel shows 
nominal stress-strain curves for the PDA (^), the PNA (a), the PAB (*), the HDA (v), 
the HNA (o) and HAB models (□). All models are evaluated at 8N = 400, b = 0.1, 
oo = 1, feeT = 1. 
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Re: Observation 1: A finite tension is required to extend an FJC from 
a corner of the unit cell to the linker, and although the vector sum of these 
tensions is zero the average tension is not zero, nor is it isotropic at the 
length scale of this unit cell. Were we to allow the boundary to relax, the 
unit cube would collapse to zero which is why the PAB model, analogous 
to what happens in classical rubber elasticity theory [20], must be stabilized 
in the ground state by a uniform hydrostatic pressure (pre-stress) to ensure 
that the cube with side oq is indeed the equilibrium volume. The value of 
this uniform pre-stress for the PAB model may be obtained by applying no 
deformation to the system (i.e., setting all Aj equal to 1 in Eq. 19). 


aps=\^/NC-^^‘' ^ 


vW 


(30) 


For the PNA and PDA it is higher, because the distributed W yield higher 
contributions to the average of (^) smaller W than from 

the larger ones rises steeply with x). We have dehned the Hookean 

spring models such that there is no pre-stress in the reference conhguration 
by setting each spring’s equilibrium length to its value in the undeformed 
state, for given linker placement. While this observation is quite obvious, 
the pre-stresses have important consequences for the non-affinity. 

Re: Observation 2: Surprisingly, the PNA model and PDA model give 
very comparable results. For FJC unit cells, the non-affinity appears negligi¬ 
ble. We verify that this is the case in Sec. This is not a universal feature: 
the equivalently dehned Hookean models do, generally, show a strong effect 
of non-affinity which, except for in uniform extension, makes the material 
softer in all deformations considered. This softening by non-affinity is a gen¬ 
eral feature, which may be attributed to the number of constraints on the 
mechanical system: the non-affine system is less constrained and therefore 
better able to accommodate strains by rotating (which costs no elastic en¬ 
ergy) rather than stretching. We propose that the suppression of non-affinity 
in the FJC models is a direct result of the pre-stress, which acts to effectively 
constrain the otherwise freely moving linker. This same effect, in extremes, 
may be seen in so-called marginal materials which may be completely floppy 
until they are stabilized by external or internal stresses |2I]. If this is the 
case, then we should be able to suppress the non-affinity by externally ap¬ 
plying a pre-stress to these unit cells as well - we demonstrate this effect in 
Sec. m 
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Re: Observation 3: The first part - that disordered FJC models are 
stiffer than the center-linker Arrnda-Boyce model, may be nnderstood by 
noting that redistribnting a hxed amonnt of chain contonr length in a nnit 
cell of hxed size necessarily yields chains shorter as well as longer than the 
average length (rather than eight chains of eqnal length). The lack of non- 
afhne displacements, noted in Observation 2, implies that the springs inside 
the nnit cnbe deform nniformly and act in parallel. The presence of non- 
afhne displacements in ensembles of springs that act in series has been noted 
before [22] . The stiffness of entropic chains scales inversely with the contonr 
length, the shorter (stiffer) chains dominate the average rigidity. Each of the 
eight springs has a spring constant 


ki 


3kBT 1 


(31) 


Comparing the two average spring constants, i.e. lJ2iki versns iJZikAB, 
after discarding constant prefactors and rewriting the nnmber of segments 
per AB chain as AAb = | Ni, we may write which holds for all partitions 
of Atot = Thns, the average and effective spring constant in the 

distribnted system is always greater than it is for the centrally placed linker. 

In Tables II and III we list the effective linear modnli extracted from onr 
simnlations. These are obtained by htting the energy graphs in Figs. |4|j^to 
a harmonic fnnction 

£(A) « ) X (A - A.J2. (32) 

where x may be a shear or a bnlk modnlus, depending on the deformation 
considered, and Agq is the A corresponding to the nnderformed state, i.e. 
Aeq = 1 for nniaxial, biaxial, pnre shear and nniform extension, and Agq = 0 
for simple shear. Table II shows that in the Hookean models, positional 
disorder (second colnmn) tends to increase all linear stiffnesses compared to 
the symmetric AB model (hrst colnmn), and that non-affinity (third colnmn) 
softens the disordered system considerably. In Appendix we demonstrate 
how the x’s for the HAB model can be determined analytically to validate onr 
method. Interestingly, the nniform extension deformation is itself symmetric, 
and absent elastic nonlinearities (as they are in the linearized regime we 
consider here) will be absorbed affinely. This explains the absence of the 
typical softening of the nniform extension modnlns between affine and non- 
affine disordered models. 
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Table 2: Linear moduli x for Hookean spring models 


Deformation 

Xhab 

Xhda 

Xhna 

Uniform extension 

36 

62.616 

61.677 

Pure shear 

0 

5.809 

0.453 

Uniaxial 

0 

4.34 

0.34 

Biaxial 

0 

17.396 

1.38199 

Simple shear 

4 

5.996 

5.245 

Table 3: Linear moduli x 

for FJC models 

Deformation 

Xpab 

Xpda 

Xpna 

Uniform extension 

38 

191.128 

188.23 

Pure shear 

48.881 

116.793 

116.095 

Uniaxial 

36.663 

87.489 

86.929 

Biaxial 

146.634 

349.889 

348.426 

Simple shear 

12.3748 

38.321 

37.605 


The differences between the FJC models and the Hookean models are in 
the amounts of pre-stress. The FJC models are unavoidably pre-stressed, 
while the Hookean models are not. These pre-stresses, we argue, serve to 
constrain otherwise free degrees of freedom within the system, and are thus 
directly related to the disappearance of the non-affinity in the FJC mod¬ 
els. We now explore this effect by applying external pre-stresses to Hookean 
systems. 

5. Pre-stress suppresses non-afRnity 

To directly measure the effect of pre-stress on non-affinity, we proceed 
now by measuring the amount of non-affinity as a function of strain, for 
differently pre-stressed systems. Several methods of quantifying the non¬ 
affinity have been proposed |22l [2S1 IMl I2SI [2S112Z| • For our systems we will 
use the differential measure proposed by Huisman et ah |25] and effectively 
measure local deviations from an affine deformation field at every strain step. 
As remarked in the previous, models with a centrally placed linker (PAB and 
HAB) will be affine by symmetry, and disordered models where we force affine 
deformations (PDA, HDA)) will, obviously, behave affinely as well. Thus we 
compare, in this section, the response of the two non-affine models, PNA 
and HNA. PNA, as we have seen, is intrinsically pre-stressed and shows little 
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(c) pure shear (d) uniform extension/compression 


Figure 6: Ensemble-averaged differential non-affine measure versus strain A for stress- 
free and pre-stressed HNA models for the applied deformation. In the main figures, the 
initial (zero-strain) lengths of individual filaments are set to 100%(*), 90%(*), 80%(B), 
70%(4), 60 %(a), 50%(<)), 40%(a), 30%(n), 20%(v), 10%(o) of the corner-to-linker dis¬ 
tance. The inset graphs show the ensemble-averaged differential non-affinity measure A^j^ 
versus strain for the PNA model, as defined in the text. 
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Figure 7: Ensemble-averaged differential non-affine measure versus strain A for stress- 
free and pre-stressed HNA models for simple shear deformation. The main figure shows 
the ensemble-averaged differential non-affinity measure A^j^ versus strain for the PNA 
model, as defined in the text. In the main figure, the initial (zero-strain) lengths of 
individual filaments are set tol00%(*), 90%(*), 80%(B), 70%(4), 60 %(a), 50%(<)), 40%(A 
), 30%(n), 20%(v), 10%(o) of the corner-linker distance. 
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non-affinity. HNA is not pre-stressed and shows non-affine softening in all 
moduli that we have measured. 

Because we consider a range of deformations, we require a slightly more 
general differential non-affinity measure than |25]. Consider, to this end, a 
linker that at some strain value A is located at X;^. Where will it be, after 
a strain increment 6X, if it moves affinely from its previous position? We 
dehne the incremental strain tensor A, as follows 

A(A + (5A) = A(A;(5A) ■ A(A). (33) 

By expanding the strain tensor around A, we obtain that 

A{\-6X) = l+(^^ ■A-'(A)^5A (34) 


For a simple shear deformation, A5 '(A;5A) is simply A5(5A), independent of 
the current state of strain. For more general deformations, however, it will 
depend on A - for instance. Ay (A; 6X) equals Ay(l -|- 5A/A). 

We construct the differential measure of non-affinity by squaring the non- 
affine deviation - he., the actual position Xa+ 5 a minus the affinely propagated 
position A(A; 6X) ■ Xa - incurred during the strain step SX, normalized to the 
magnitude of the strain step: 


A 




Xa+5a-A(A;5A)-Xa 


(35) 


This quantity is plotted against the strain A, for all deformation modes con¬ 
sidered, for the two non-affine models (PNA and HNA) in Figs. [6]j^ For the 
HNA models, we apply an external pre-stress by reducing the rest lengths of 
all springs by a uniform factor. The resulting network will generally require 
a second relaxation towards a new ground state. We strain the systems from 
this new, equilibrated, conhguration to obtain the measures displayed here. 
All HNA graphs clearly demonstrate that indeed, also in the Hookean models 
pre-stress results in diminished non-affinity. What these graphs also show, is 
that although the non-affinity is small for the FJC models, it is not zero. No¬ 
tably, even the non-prestressed Hookean systems show a continued downward 
trend in differential non-affinity at hnite strain, whereas the FJC systems in 
all cases show a rise over the regime considered. The differential measure is 
expected to tend to zero in the large strain limit. At the highest forces, the 
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linker is completely constrained and no longer free nor able to move other 
than affinely. This high strain limit, however, may not be attained in the 
FJC models as these are intrinsically hnitely extensible: it is impossible to 
extend them beyond a contour length of Nb. We postulate that in our sys¬ 
tems, the chains are too short to attain the affine limiting behavior. Our 
Hookean systems show a decreasing non-affinity, while the FJC models show 
a small but rising non-affinity. While there are very few studies that measure 
the strain-dependence of the non-affinity, both behaviors have actually been 
reported. The increase of non-affininty is seen in simulations on Worm Like 
Chain systems, which also have pronounced non-linear elasticity but not, 
necessarily, pre-stresses. Fibrin gels, in the only experimental measurement 
that we are aware of, display a pronounced peak in non-affinity before en¬ 
tering a prolonged regime of diminishing non-affine deformations [28]. We 
see, that depending on the deformation mode considered the amount of non¬ 
affinity varies greatly - it appears that some deformations are more conducive 
to non-affinity than others. 

6. Summary and Conclusions 

We set out to assess what the effect of spatial and elastic disorder, single¬ 
chain elasticity and pre-stress is on the mechanical response of simple unit 
cell network models for rubber-like elasticity. Models based on the Freely 
Jointed Chain are invariably pre-stressed. Addition of a uniform pressure 
stabilizes a state, but does not remove pre-stress. FJC models are intrin¬ 
sically less non-affine than Hookean spring models, even when these are, 
otherwise, elastically equivalent. The effects of externally applied pre-stress 
in Hookean systems are similar, as they too suppress the amount of non¬ 
affinity, but closer inspection reveals that the manner in which non-affinity 
depends on strain is still quite different. We speculate that this is due to the 
nonlinear force extension relation of the FJC. 

In reality, many materials are indeed pre-stressed. For biopolymer gels 
and tissues, in particular, the presence of pre-stresses is well-documented. 
Previous modeling [22] assumed complete affinity, yet was able to capture 
the mechanical response of a variety of biopolymer gels reasonably well. A 
possible explanation for the success of such affine models, we propose, may be 
that non-linearly elastic polymer models, even in the most basic incarnations 
considered here, display a greatly suppressed non affinity and are, by many 
accounts, very well approximated by fully affine models. 
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This finding also snggests one should proceed with caution in replacing 
nonlinear constituents with linear springs. If one removes the pre-stress, this 
may result in an unrealistically high degree of non-affinity and, accordingly, 
lower network stiffness. Even if one linearizes the springs retaining an ap¬ 
propriate amount of pre-stress, the non-affine response may still evolve very 
differently with strain due to the nonlinearities. 

In those experiments that are able to map out the non-affinity of the 
crosslinkers (or of otherwise fixed fiducial markers inside the network) [221 
l2n [28l [30l 1^ . it should be feasible to investigate directly the correlation 
between pre-stress and non-affinity. We predict these two to be inversely 
related. 

In upcoming work, we explore the effects of polymer persistence on the 
effects reported here. The associated Worm-Like Chain model features a 
force-extension relation that has both a finite rest length, as well as a steeply 
nonlinear response to stretching. As such, semifiexible chains represent an 
interesting middle ground between the Hookean and FJC cases considered 
here. It is entirely unclear which effect will be dominant in their mechanical 
response: the lack of ground state pre-stresses promoting non-affinity, or the 
nonlinear force-extension that suppresses it. 
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Appendix A. Numerical computation of moduli and exact expres¬ 
sions for the HAB model 

As detailed in the text, we compute moduli by determining, numerically, 
the curvature of the energy vs. strain graph. In one case considered here, we 
can compute this exactly. For the HAB model, we have eight chains that each 
contribute identically to the strain energy. One of them, in the undeformed 
state, has as its end-to-end vector v = {ao/2){x + y + z), and its length 
is \/3ao/2 [13]. Acting on this polymer with, for instance, the deformation 
tensor for simple shear, it is transformed affinely to v' = /\s{X)-v. The length 
of this deformed vector is 


— —cio'\/J- 2A J- 3 . 


(A.l) 



Writing the deformation energy of each chains as (l/2)/csp5£^, multiplying 
by 8 for the number of chains in the unit cell, and dividing by the unit cell 
volume we arrive at the strain energy density for the HAB in simple shear 


^5(A) = ^ fVA2 + 2A + 3-^)" . 

ao V ) 

To extract the modulus, we expand to second order in A: 
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From which we read off the modulus in simple shear 

2L 


Xs = 




3ac 
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(A.4) 


For the values we have chosen in the simulations (W = 50,6 = = 

l,ao = 1, the spring constant computed with Eq. 11 is fcgp = 6 and thus 
X = 4. This was conhrmed in the simulations. All other moduli may be 
similarly computed - it turns out that the only other nonzero one is the bulk 
modulus probed in uniform extension; xv = (Sfcsp/oo) = 36. 
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